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C ! Abstract 
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Quantum-mechanical "PT-symmetric theories associated with complex cubic 
potentials such as V = x 2 + y 2 + igxy 2 and V = x 2 + y 2 + z 2 + igxyz, where 
g is a real parameter, are investigated. These theories appear to possess 
real, positive spectra. Low-lying energy levels are calculated to very high 
order in perturbation theory. The large-order behavior of the perturbation 

o 

■ coefficients is determined using multidimensional WKB tunneling techniques. 

This approach is also applied to the complex Henon-Heiles potential V = 
x 2 + y 2 + igixy 2 — 5X 3 ). 
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^ . In this Letter we examine complex PT-symmetric cubic Hamiltonians such as 
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H (2) = p 2 x + p 2 y + x 2 + y 2 + igxy 2 
H {3) =p 2 x +p 2 y +p 2 z + x 2 + y 2 + z 2 + igxyz, (1) 

where g is a real parameter [|lj. The superscript on H indicates the number of degrees 
of freedom; these Hamiltonians are several-degree-of-freedom generalizations of the one- 
dimensional complex cubic Hamiltonian = p 2 + x 2 + igx 3 , which has recently been 
studied in great detail by many authors fl^l2|. The non-Hermitian Hamiltonian is 



interesting because its spectrum is entirely real and positive. The reality of the spectrum is 
apparently due to the VT invariance of the Hamiltonian. 



While many different one-degree-of-freedom examples of non-Hermitian PT-symmetric 
quantum systems have been studied, no multidimensional complex PT-symmetric coupled- 
oscillator systems have been examined. The purpose of this Letter is to show that (i) the 
property of real, positive spectra persists even for quantum systems having several degrees 
of freedom, and (ii) these theories have many other properties in common with theories 
described by conventional Hermitian Hamiltonians. 

Direct numerical evidence for the reality and positivity of the spectrum of can 
be found by performing a Runge-Kutta integration of the associated complex Schrodinger 
equation 0] . Alternatively, the large-energy eigenvalues of the spectrum can be calculated 
with great accuracy by using conventional WKB techniques [[HJ. A strong argument for 
the reality and positivity of the spectrum can be obtained by calculating the spectral zeta 
function Z(l) (the sum of the inverses of the eigenvalues). For the Hamiltonian 1 H = p 2 +ix 3 
this was done by Mezincescu || and Bender and Wang [10]. The exact result for Z(l) is 



4sin 2 (7r/5)r 2 (l/5) 

Z[l} ~ 5*V5r(3/5) • {2) 

Using the numerical values of the first few eigenvalues and the WKB formula for the high 
eigenvalues, one can conclude that any complex eigenvalues must be larger in magnitude 
than about 10 18 . Some rigorous results regarding the reality of the eigenvalues of H^ l > have 



been obtained by Shin [JTTJ] , who showed that the entire spectrum must lie in a narrow wedge 
containing the positive-real axis. Other results have been obtained by Delabaere et al ||,[7j . 

Let us now return to the Hamiltonians in ([!]). The Schrodinger equations associated 
with H^ 2 \ and are 

— i/)"(x) + (x 2 + igx 3 )?pi(x) = Eifti(x), 



1 This Hamiltonian, the massless case of has a positive discrete spectrum. It is not known if 
the massless versions of and have discrete spectra. Indeed, even for the massless coupled 
anharmonic oscillator potential V = x 2 y 2 , it is also not known if the spectrum is discrete. 
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-V 2 ip 2 (x, y) + {x 2 + y 2 + igxy 2 )ip 2 (x, y) = Eip 2 (x, y), 
-V 2 ip 3 (x, y, z) + (x 2 + y 2 + z 2 + igxyz)ip 3 (x, y, z) = Etjj 3 (x, y, z). (3) 

We have solved the Schrodinger equations (|3]) for the eigenvalues in several ways. One 
technique is to diagonalize each Hamiltonian in a set of multidimensional harmonic oscillator 
basis states. This procedure immediately reveals that the energy levels are real. 

A more precise calculation of the energies of the complex "PT-symmetric Hamiltonians in 
(H) is performed using high-order Rayleigh-Schrodinger perturbation theory. This technique 



was used in Refs. || and to obtain the perturbation series for the ground-state energy 
of In these references it was found that the Rayleigh-Schrodinger perturbation series 

is Borel summable and that Pade summation is in excellent agreement with the real energy 
spectrum. Furthermore, Pade analysis provides strong numerical evidence that the once- 
subtracted ground-state energy considered as a function of g 2 is a Stieltjes function. 

The Rayleigh-Schrodinger perturbation series for the ground-state energies Eq~\ Eq 2 \ 
and -Eg 3 " 1 of the Hamiltonians H^ 2 \ and have the asymptotic form 

n) 11 „ 465 , 39709 , 19250805 „ 2944491879 in 

EX ' ~ H q 2 g 4 H o 6 q 8 H g 10 H , 

16 y 256 y 4096 y 262144 y 4194304 y 

r 2) 5 2 223 . 114407 „ 346266143 R 2360833242959 in 

Ek ' ~ 2 + — g 2 o 4 + o 6 g 8 + o 10 + • • • , 

48 6912 y 4976640 14332723200 72236924928000 

(3) - , 1 2 7 4 5069 6 2441189 8 

Ek y ~ 3 H q q H q q 

48 y 4608 y 19906560 y 38220595200 y 

8034211571 1(1 

H 9 H (4) 

385263599616000 y v ; 

in the limit g — > 0. The (9, 9) Pade was constructed from the once-subtracted form of these 
series and the results are plotted in Fig. |I[ The first two excited states of are plotted 
in Fig. |[ Note that the degenerate unperturbed energy level at E = 4 splits into two levels, 
each of which is greater than four. 

The perturbation coefficients in these series are derived from recursion relations like 



those first derived for the anharmonic oscillator []TJJ. These recursion relations are obtained 
directly from the Schrodinger equations (^). We substitute ip\(x) = e~ x2 ^ 2 (pi(x), ^{x^y) = 
e~( x2+y2 ^ 2 (j)2(x,y) and ip 3 (x,y,z) = e~ ( - x2+y2+z2 ^ 2 (j) 3 (x,y, z), where 0i,2,3 are formal power 



series in g; the coefficients of g n are polynomials P n of degree 3n in the variables x, y, z. 
For example, for the case of r/> 3 , P n = Y^j^,i=Q a n,j,k,ix^y k z l , and the coefficients a n ,j,k,i satisfy 

1 



(^n,j,k,l 



n-1 



0"n-l,j-l,k-l,l-l — 2 ^ Ctn-pj,fc,/(flp,2,0,0 + a p,0,2,0 + %>,0,0,2) 
p=l 



2(j + k + 1) 

+(j + I) (J + 2)a n> j + 2,k,i + (k + l)(k + 2)a nj - fc +2,« + (/ + !)(/ + 2)a„j ife ,j +2 



• (5) 



Once the coefficients a n ,j,k,i are known, we can construct the coefficient of g 2n in the expansion 
for the ground-state energy Eq in (B) according to 



E^ 3) ~ 3 + 2 ^(a 2ni2 ,o,o + a2n,o,2,o + a2n,o,o,2)(-fl ,2 ) n (# -»• 0). 



(6) 



n=l 



We are particularly interested in the large-order behavior of the coefficients in the per- 
turbation expansion because this behavior suggests that the series is Borel summable and 
reveals the analytic structure of the energy level as a function of complex g 2 . In Ref. it 



is shown that the large-n behavior of the coefficient of g 2n in the expansion of Eq 1 ^ is 



4 /l^\ n+1 / 2 / 1 



l-0(- 

n 



in — > oo . 



(7) 



Therefore, although divergent, the series for Eq 1 ^ is Borel summable ||13|| . Observe that if 
the factor of i were absent from the Hamiltonian then the perturbation coefficients 

would not alternate in sign and the perturbation series would not be Borel summable. 
A major result reported here is the large-order behavior of the coefficients of g 2n in the 



series for E, 



(2) 
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xn+1 



72V2 /j.yVn + i 
TrVcoshdTrv 7 ^) U87 V 2 



1-0 



n 



in — > oo 



(8) 



and coefficients of g 2n in the series for Ef 



-1) 



n+l 



1152^ 



5 r +1/2 r ^ + i 



1-0 



/? 



(n — > oo). 



(9) 



v^F cosh(|7r\/23) \72j 

We have verified these results to extremely high precision by performing a Richardson ex- 
trapolation |]I3| of the perturbation coefficients in (||) divided by these behaviors. 



We derive these results by adapting the multidimensional WKB tunneling techniques in 
Ref. [[DJ. We observe that if g is replaced by ig in and then we obtain potentials 
for which the probability current in a Gaussian ground state leaks out to infinity. The 
probability flows outward along most-probable escape paths (MPEPs). 

To determine the MPEPs for we rewrite the potential in polar coordinates: 

V(x, y) = x 2 + y 2 — gxy 2 = r 2 — gr 3 cos 9 sin 2 9, (10) 

where x = r cos 9, y = r sin 9. Letting a = g cos 9 sin 2 9, we calculate V r = 2r — 3ar 2 . Then, 
setting V r — gives the critical radius r = and at this radius the potential has the 
value -§fOt~ 2 . Thus, V achieves its minimum when sin# = ±y§- Therefore, the effective 
radial potential is V(r) = r 2 — i^r 3 . We conclude that there are two straight-line MPEPs 
symmetrically placed above and below the positive-x axis (see Fig. |3|). 

Geometrical optics (ray tracing) is sufficient to reproduce the gamma-function and ex- 
ponential behaviors in (pi). We simply evaluate the approximate WKB integral 



3\/3 
29 



-2 / *" dr a r 2 ^=r 3 , 



o V 3^ 

where we have neglected the constant term in the limit of small g. We evaluate the resulting 
beta-function integral to get the leading exponent in the tunneling rate: / = —^g~ 2 . There 
is a standard dispersion-integral procedure |6|JI^| that expresses the large-order behavior as 



the nth inverse moment of the tunneling rate. This procedure gives the behavior in @ 
apart from an overall multiplicative constant. This constant can only be determined by 
performing a physical- optics calculation of the tunneling rate. 

For this physical-optics calculation we must determine the flux of probability through a 
tube centered about the MPEP. We introduce a rotated coordinate system by 

1 y/2 + 1 , y/2 

x = —j=r pi, y = —=t H — P r, 

so that r measures the distance along the MPEP and t is the coordinate transverse to the 
MPEP. In terms of these variables, the Schrodinger equation @ for ip2 reads 



- V 2 ^M) + 



r 2 + t 2 + -^(-2r 3 + 3rt 2 - v 7 ^ 3 ) - 2 



^ 2 (r,t)=0, (11) 



3v/3 

where we have replaced g by ig. We may drop the t 3 term because gt 3 «t 2 for small g. 

Next, we separate the radial dependence from the transverse dependence by writing 
^2 (?">£) = W(r)<p(r, t). The function W(r), which expresses the radial dependence, satisfies 
the differential equation W"{r) = (r 2 — ^/g 7 " 3 — l) W(r). The WKB approximation to the 
decaying solution to this equation is 

^ ('- 2 - ^ - 

where the numerical factors are included in anticipation of asymptotic matching. The equa- 
tion for <f>(r,t) is -2^<p r - <p rr - <p tt + (t 2 + -^rt 2 - 1)0 = 0. Note that W r /W r for 



small g. The change of variable v = — yields a parabolic equation for </>: 

(v 2 - 1)0, - <f> a + [t 2 + 3t 2 (l - t; 2 )/2 - 1]0 = 0, (13) 

where we neglect the small term of order g 2 (j) vv . The solution to this equation has the form 
of a Gaussian that expresses the thickness of the stream of probability current that flows 
outward along the MPEP: 

<f>(v,t) = A{v)e- t2f{v)/2 , (14) 
where the function f(v) satisfies the Riccati equation 

(l-v 2 )f(v)-2f 2 (v) + 5-3v 2 = (15) 
and A(v ) satisfies the transport equation 

(1 - v 2 )A'{v) - f{v)A{v) + A{v) = 0. (16) 

To solve the Riccati equation fli~5|) we substitute /(f) = — — v 2 )h'{y )/h{v ) and convert 
it to the second-order linear equation 

(1 - v 2 )h"{v ) - 2vh'{v) + (-G - K v ) = 0, (17) 
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which we recognize as the Legendre differential equation [16 



To find the initial conditions on f(v) and A(v), we match <fr(v,t) to the wave-function 
solution to the Schrodinger equation flTT|) in the inner region where r and t are of order 1. In 
this region the wave function is a Gaussian: (p(v,t) = e~^ r2+t2 ^ 2 . By construction, for small 
r, W(r) ~ e -r2 / 2 . Thus, we obtain the initial conditions A(l) = 1 and /(l) = 1. Hence, the 
solution to (|17D is the Legendre function h{v) = P~ 2 (v), where z/(z/ + 1) = —6. 



Our objective now is to find the flux of probability at the distant turning point ||15|| . At 
this turning point the radial component of the probability current is 



1 - *£j?* m *w °* (- 2 Ir ds ft - ih 8 ' 3 - j • (18) 

where we have included a factor of 2 because there are two channels. We integrate in the 



transverse direction to get the total flux of probability current /f^, dt e 1 ^ = ^Jj^yj and 
we evaluate the integral in the exponent to obtain 

2 l V da ft ~ iTT 3 " 1 ~ W ~ " l(12v/I) + ln(9) " 5 (9 ^ 0) ' 

Thus, the total outward flux of probability is ^U 2 (0) x f^e' 18/ ^ g2) . Substituting this re- 



suit into the dispersion integral, we obtain the following formula for the large-order behavior 
of the coefficients in the perturbation series for the ground-state energy: 



ir 12 ^flV +1/ V n + l - 



V18 



1 - O 



n 



In 



oo). (19) 



It remains to find the numbers /(0) and A(0). Using the hypergeometric-function rep- 
resentation for the Legendre function [[TjJ, we obtain 

. . . T(u/2)T(2 + u/2) 

f = - tan 7TZ//2 — , v ' 1 , \ , ' — — 

JK J v 1 ' r(-l/2 + i//2)r(3/2 + i//2) 

and 



A(0) = 7T- 1 / 4 J2T{2 + z//2)r(3/2 - v/2). 



Thus, A 2 (0)/Jf(0) = v /24vr/cosh(7rv / 23/2) and we have derived the result in 
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To obtain the formula (H) we follow the same procedure. In this case there are four 
radial MPEPs and there are two transverse variables. We begin by introducing a change of 
coordinates in the Schrodinger equation (Q) for ip 3 (x,y, z): 

111 111 12 
x = —;=r H -;=s i=t, y = —=r H y=s H =t, z = — 

The remainder of the calculation is identical to that summarized above for ip 2 {x,y)- 
We conclude by noting that the complex Henon-Heiles potential 

V HH = x 2 + y 2 + ig{xy 2 - \x 3 ) (20) 

has a Borel summable Ray leigh- Schrodinger perturbation series for the ground-state energy, 
with the leading growth of the coefficients given by 



c(-i) n+1 |r) r(n+ i 



24/ V 2 



1 + 



(n-oo), (21) 



where C is a constant. A simple geometric-optics calculation (with g replaced by ig) con- 
firms this leading growth rate. There are three MPEPs, a pair of MPEPs similar to those 
encountered in our analysis of (now at angles ±| from the positive-x axis), and a third 
MPEP along the negative-x axis, which is like that for (see Fig. ^). Remarkably, these 
two different types of MPEPs produce exactly the same leading contribution to fl2jj) . This 
fact depends crucially on having the appropriate combinatorial factors in (f2~0"D. 

Finally, we remark that the quantum field theoretic generalizations of the Hamiltonians 
studied here, particularly H^, may be viewed as theories of scalar electrodynamics. It 
would be interesting to study such issues as bound states and Schwinger-Dyson equations 
in such theories. 
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FIG. 1. Ground-state energies of the Hamiltonians H^> (solid line), (long-dashed line), 

and (short-dashed line), as functions of the coupling constant g. Note that the energy lev- 
els are real and positive. The graphs were obtained from the (9, 9) Pade constructed from the 
once-subtracted perturbation series for these energy levels. 
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FIGURE 2 



energy 




FIG. 2. First two excited energies of the Hamiltonian . The solid and dashed lines represent 
the levels whose unperturbed states are xe 

-(xW)/2 and ye -(x 2 +y 2 )/2_ The un p er t ur bed energies 
split into two distinct levels, both of which lie above the unperturbed level at E = 4. The graphs 
were constructed from the (9, 9) Pade of the perturbation expansion in powers of g 2 . 
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FIGURE 3 




-3-2-10 1 2 3 



FIG. 3. Contour plot of the potential x 2 +y 2 —xy 2 . A Gaussian probability distribution localized 
at the origin gradually leaks out to infinity preferentially along two channels in the right-half plane. 
These channels are called most probable escape paths (MPEPs). 
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FIGURE 4 
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FIG. 4. Contour plot of the Henon-Heiles potential x 2 + y 2 — xy 2 + ^x 3 . A Gaussian probability 
distribution localized at the origin tunnels out to infinity preferentially along three MPEPs, two 
in the right-half plane and one along the negative-real axis. Remarkably, the contribution from all 
three MPEPs to the large-order behavior of perturbation theory is of the same magnitude. 



14 



